Reference site numbering

So far all the preceding examples have involved analysis of proteins with sequential integer site numbers (e.g., 1, 2, 3, …). This sequential integer numbering can also start at sites other than one (e.g., 331, 332, 333, …).

However, often proteins are numbered in relation to alignment to sequential integer numbering of some reference homolog. For instance, SARS-CoV-2 spike sequences are usually numbered in relation to the Wuhan-Hu-1 reference. But if the protein has accumulated insertions or deletions relative to the reference, there may be missing sites (gaps) or insertions (typically numbered as 214, 214a, 214b, etc).

This notebook shows how to perform an analysis and analyze the results with non-sequential reference site numbering. The advantage of doing this is that the results can directly be visualized/analyzed using the conventional site numbering scheme.

It does this by analyzing deep mutational scanning of the SARS-CoV-2 Omicron BA.1 spike, which has several indels relative to the Wuhan-Hu-1 numbering reference.

Fit and visualize a model with non-integer site numbering

First, import Python modules:

import polyclonal

import pandas as pd

Now read the data to fit:

# read data w `na_filter=None` so empty aa_substitutions read as such rather than NA
data_to_fit = pd.read_csv(
    "Lib-2_2022-06-22_thaw-1_LyCoV-1404_1_prob_escape.csv",
    na_filter=None,
)

data_to_fit.head()

	antibody_concentration	barcode	aa_substitutions_sequential	aa_substitutions_reference	n_aa_substitutions	prob_escape	prob_escape_uncensored	no-antibody_count
0	0.654	AAGAGCTTACTCTCGA	S443P S1167P	S446P S1170P	2	0.9004	0.9004	2890
1	0.654	ATAAGATAGATTTAGG	K441T A519S	K444T A522S	2	0.8006	0.8006	2217
2	0.654	GATCGAGTGTGTAGCA	F152T K441E	F157T K444E	2	0.5797	0.5797	2966
3	0.654	CCAAACGGTATGATGA	Q14H V67P D212H N445T I1224M	Q14H V67P D215H N448T I1227M	5	0.4125	0.4125	4131
4	0.654	TTACTGTGCAACCCAA	F163L N445S	F168L N448S	2	0.3600	0.3600	4275

Notice how these data have sites numbered in two different schemes:

1. aa_substitutions_sequential: sequential integer numbering (1, 2, 3, …) of the protein used in the experiment. You could just analyze the mutations this way, but then the results will not be in standard reference numbering scheme.

2. aa_substitutions_reference: the referencey based site numbering, which skips sites with indels and has some non-numeric sites where there are insertions (eg, 214a), as for example in the variant shown below:

data_to_fit.query("aa_substitutions_reference.str.contains('a')").head(n=1)

	antibody_concentration	barcode	aa_substitutions_sequential	aa_substitutions_reference	n_aa_substitutions	prob_escape	prob_escape_uncensored	no-antibody_count
32	0.654	CTAGATAAATCCCTGC	E209A K441N	E214aA K444N	2	0.3965	0.3965	2214

Here we will use the reference based amino-acid substitutions, so create a column with the name used by Polyclonal (the “aa_substitutions” column) that uses this numbering scheme:

data_to_fit["aa_substitutions"] = data_to_fit["aa_substitutions_reference"]

Importantly, in order to use reference based numbering that is not sequential integer, you also have to provide a list of the sites in order. The reason is that otherwise it’s not possible for Polyclonal to figure out for instance if sites are just missing from the data are are actually deletions.

Here we read a data frame that maps sequential to reference site numbering, and then use that to extract the list of reference sites:

site_numbering_map = pd.read_csv("BA.1_site_numbering_map.csv")
display(site_numbering_map.head())

print("Note how some reference sites differ from sequential ones due to indels:")
display(
    site_numbering_map[
        site_numbering_map["sequential_site"].astype(str)
        != site_numbering_map["reference_site"]
    ].head()
)

sites = site_numbering_map["reference_site"].tolist()

	sequential_site	reference_site
0	0	-1
1	1	1
2	2	2
3	3	3
4	4	4

Note how some reference sites differ from sequential ones due to indels:

	sequential_site	reference_site
0	0	-1
69	69	71
70	70	72
71	71	73
72	72	74

Now initialize and fit the Polyclonal model, but pass the sites argument so we can use these non-sequential-integer reference sites. (If you don’t pass the sites argument, sites are assumed to be sequential integer):

model = polyclonal.Polyclonal(
    # `polyclonal` expects the concentration column to be named "concentration"
    data_to_fit=data_to_fit.rename(columns={"antibody_concentration": "concentration"}),
    n_epitopes=1,
    alphabet=polyclonal.AAS_WITHSTOP_WITHGAP,
    sites=sites,
)

Note how the model has its sequential_integer_sites attribute set to False:

assert set(model.sites) == set(sites)

assert model.sequential_integer_sites is False

Now fit the model:

# NBVAL_IGNORE_OUTPUT

_ = model.fit(logfreq=200, reg_escape_weight=0.1)

# First fitting site-level model.
# Starting optimization of 1248 parameters at Tue Aug  9 06:16:17 2022.
         step     time_sec         loss     fit_loss   reg_escape   reg_spread reg_activity
            0     0.052611        37315        37314            0            0      0.90499
          188       11.514       4865.6       4822.7       39.237            0       3.6724
# Successfully finished at Tue Aug  9 06:16:29 2022.
# Starting optimization of 8450 parameters at Tue Aug  9 06:16:29 2022.
         step     time_sec         loss     fit_loss   reg_escape   reg_spread reg_activity
            0     0.077885       7299.1       7031.4       264.09   8.7534e-31       3.6724
          200       17.083       6821.7       6715.9       80.416       21.021       4.3081
          299       25.066       6821.1       6714.3       81.169       21.295       4.3101
# Successfully finished at Tue Aug  9 06:16:54 2022.

Look at output. First we look at the mut_escape_df. The site entries are str:

# NBVAL_IGNORE_OUTPUT

assert all(model.mut_escape_df["site"].astype(str) == model.mut_escape_df["site"])
assert set(model.mut_escape_df["site"]).issubset(sites)

model.mut_escape_df.head().round(2)

	epitope	site	wildtype	mutant	mutation	escape	times_seen
0	1	1	M	I	M1I	-0.00	1
1	1	1	M	T	M1T	-0.01	2
2	1	2	F	I	F2I	0.00	2
3	1	2	F	L	F2L	0.01	14
4	1	2	F	S	F2S	-0.01	14

The same is true for mut_escape_site_summary_df:

# NBVAL_IGNORE_OUTPUT

assert all(
    model.mut_escape_site_summary_df()["site"].astype(str)
    == model.mut_escape_site_summary_df()["site"]
)
assert set(model.mut_escape_site_summary_df()["site"]).issubset(sites)

model.mut_escape_site_summary_df().head().round(2)

	epitope	site	wildtype	mean	total positive	max	min	total negative	n mutations
0	1	1	M	-0.00	0.00	-0.00	-0.01	-0.01	2
1	1	2	F	0.01	0.03	0.02	-0.01	-0.01	4
2	1	3	V	0.06	0.64	0.56	-0.12	-0.14	8
3	1	4	F	0.22	1.12	1.02	-0.01	-0.01	5
4	1	5	L	0.02	2.08	0.94	-0.89	-1.78	20

Now plot the mutation-escape values:

# NBVAL_SKIP
# NBVAL_IGNORE_OUTPUT
# we ignore this cell with `nbval` because for some reason it hangs on Travis

model.mut_escape_heatmap(init_min_times_seen=3)

And plot the site-escape values:

# NBVAL_IGNORE_OUTPUT

model.mut_escape_lineplot(min_times_seen=3)

Confirm same results for reference and sequential integer number

To demonstrate how the results are the same regardless of which numbering scheme is used for the fitting, we also fit a model with the sequentially numbered variants. This section of the notebook can almost be considered a test rather than an example.

# NBVAL_IGNORE_OUTPUT

model_sequential = polyclonal.Polyclonal(
    data_to_fit=(
        data_to_fit.drop(columns="aa_substitutions").rename(
            columns={
                "antibody_concentration": "concentration",
                "aa_substitutions_sequential": "aa_substitutions",
            }
        )
    ),
    n_epitopes=1,
    alphabet=polyclonal.AAS_WITHSTOP_WITHGAP,
)

assert model_sequential.sequential_integer_sites is True

_ = model_sequential.fit(logfreq=200, reg_escape_weight=0.1)

# First fitting site-level model.
# Starting optimization of 1248 parameters at Tue Aug  9 06:17:01 2022.
         step     time_sec         loss     fit_loss   reg_escape   reg_spread reg_activity
            0     0.050311        37315        37314            0            0      0.90499
          188       11.138       4865.6       4822.6       39.396            0       3.6721
# Successfully finished at Tue Aug  9 06:17:12 2022.
# Starting optimization of 8450 parameters at Tue Aug  9 06:17:12 2022.
         step     time_sec         loss     fit_loss   reg_escape   reg_spread reg_activity
            0      0.06758         7300       7031.3       265.04   8.9894e-31       3.6721
          200       15.792       6821.6       6715.8       80.443       21.032       4.3093
          245       19.444       6821.4       6715.2       80.692       21.171       4.3108
# Successfully finished at Tue Aug  9 06:17:32 2022.

Now make sure the fitting gives nearly the same result regardless of which site numbering is used. First check the activity values:

pd.testing.assert_frame_equal(
    model_sequential.activity_wt_df,
    model.activity_wt_df,
    atol=0.01,
)

Now compare the mutation-escape values. In order to do this, we have to re-number the sequential values to reference numbering:

min_times_seen = 50

mut_escape = model.mut_escape_df.drop(columns="mutation")

mut_escape_sequential = model_sequential.mut_escape_df.assign(
    site=lambda x: x["site"].map(
        site_numbering_map.set_index("sequential_site")["reference_site"].to_dict()
    )
).drop(columns="mutation")

# have to use fairly big atol to test this, so also do correlations
pd.testing.assert_frame_equal(
    mut_escape,
    mut_escape_sequential,
    atol=1.5,
)
assert 0.99 < mut_escape["escape"].corr(mut_escape_sequential["escape"])

Test of averaging

Just as a toy test, average the reference-site model with itself and make sure it still displays reference-based site numbers:

avg_model = polyclonal.PolyclonalAverage(
    pd.DataFrame({"number": [1, 2], "model": [model, model]}),
)

# NBVAL_IGNORE_OUTPUT

avg_model.mut_escape_heatmap(init_min_times_seen=3)

# NBVAL_IGNORE_OUTPUT

avg_model.mut_escape_lineplot()

You cannot average models that don’t have the same site-numbering setup in terms of sequential or non-sequential:

# NBVAL_RAISES_EXCEPTION

# this will give an error
polyclonal.PolyclonalAverage(
    pd.DataFrame({"number": [1, 2], "model": [model_sequential, model]}),
)

---------------------------------------------------------------------------
PolyclonalHarmonizeError                  Traceback (most recent call last)
/tmp/ipykernel_62444/345686643.py in <module>
      2
      3 # this will give an error
----> 4 polyclonal.PolyclonalAverage(
      5     pd.DataFrame({"number": [1, 2], "model": [model_sequential, model]}),
      6 )

~/polyclonal/polyclonal/polyclonal_collection.py in __init__(self, models_df, harmonize_to, default_avg_to_plot)
    920             harmonize_to = models_df.iloc[0]["model"]
    921
--> 922         models_df["model"] = [
    923             m.epitope_harmonized_model(harmonize_to)[0] for m in models_df["model"]
    924         ]

~/polyclonal/polyclonal/polyclonal_collection.py in <listcomp>(.0)
    921
    922         models_df["model"] = [
--> 923             m.epitope_harmonized_model(harmonize_to)[0] for m in models_df["model"]
    924         ]
    925

~/polyclonal/polyclonal/polyclonal.py in epitope_harmonized_model(self, ref_poly)
   2160
   2161         if self.sequential_integer_sites != ref_poly.sequential_integer_sites:
-> 2162             raise PolyclonalHarmonizeError(
   2163                 "Models being harmonized don't have same `sequential_integer_sites`"
   2164             )

PolyclonalHarmonizeError: Models being harmonized don't have same `sequential_integer_sites`